SEMSummary <- function(formula, data) {
  vars <- attr(terms(formula, data = data), "variables")
  vnames <- as.character(vars)[-1L]
  if (length(vnames) < 2) stop("You must specify at least 2 variables to use this function")
  env <- environment(formula)

  if (!is.data.frame(data)) {
    data <- as.data.frame(data)
  }

  X <- eval(vars, data, env)
  names(X) <- vnames
  X <- as.data.frame(X)

  rm(data)
  gc()

  mu <- colMeans(X, na.rm = TRUE)
  stdev <- sapply(X, sd, na.rm = TRUE)
  Sigma <- cov(X, use = "pairwise.complete.obs")
  sSigma <- cor(X, use = "pairwise.complete.obs")

  n <- nrow(X)
  L <- is.na(X)
  nmiss <- colSums(L)
  i <- which(upper.tri(Sigma), arr.ind = TRUE)
  pairmiss <- apply(i, 1L, function(j) {
    sum(L[, j[1]] | L[, j[2]])
  })
  coverage <- matrix(NA, nrow = ncol(X), ncol = ncol(X))
  diag(coverage) <- (n - nmiss)/n
  coverage[i] <- (n - pairmiss)/n
  coverage[i[, c(2, 1)]] <- (n - pairmiss)/n
  dimnames(coverage) <- dimnames(Sigma)


  names(nmiss) <- names(mu) <- names(stdev) <- names(X)

  output <- list(names = vnames, n = n, nmissing = nmiss, mu = mu, stdev = stdev,
    Sigma = Sigma, sSigma = sSigma, coverage = coverage)
  class(output) <- "SEMSummary"

  return(output)
}

